

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>

#define alpha 1
#define beta 1
double A[N][N +20];
double B[N][N +20];
double x[N];
double u1[N];
double u2[N];
double v1[N];
double v2[N];
double w[N];
double y[N];
double z[N];

void init_arrays()
{
    int i, j;
    for (i=0; i<N; i++) {
        u1[i] = i;
        u2[i] = (i+1)/N/2.0;
        v1[i] = (i+1)/N/4.0;
        v2[i] = (i+1)/N/6.0;
        y[i] = (i+1)/N/8.0;
        z[i] = (i+1)/N/9.0;
        x[i] = 0.0;
        w[i] = 0.0;
        for (j=0; j<N; j++) {
            A[i][j] = ((double) i*j)/N;
        }
    }
}

double rtclock()
{
  struct timezone tzp;
  struct timeval tp;
  int stat;
  gettimeofday (&tp, &tzp);
  return (tp.tv_sec + tp.tv_usec*1.0e-6);
}

int main()
{
  init_arrays();

  double annot_t_start=0, annot_t_end=0, annot_t_total=0;
  int annot_i;

  for (annot_i=0; annot_i<REPS; annot_i++)
  {
    annot_t_start = rtclock();



/*@ begin PerfTuning (
  def build
  {
  arg command = 'icc';
  arg options = '-fast -openmp -I/usr/local/icc/include -lm';
  }

  def performance_counter
  {
  arg method = 'basic timer';
  arg repetitions = 4;
  }  
  
  def performance_params
  {  
    param T1_1[] = [1];
    param T1_2[] = [1];
    param T2_1[] = [1];
    param T2_2[] = [1];

    param U1[] = [1];
    param U2[] = [4];
    param U3[] = [4];
    param U4[] = [1];
    param U5[] = [4];

    param PERM1[] = [
#      [0,1],
      [1,0],
    ];
    param PERM2[] = [
      ['i','j'],
#      ['j','i'],
    ];

    param PAR[] = [False];
    param SCREP1[] = [True];
    param IVEC1[] = [True];
    param SCREP2[] = [True];
    param IVEC2[] = [True];
  }
  
  def search
  {
  arg algorithm = 'Exhaustive';
# arg algorithm = 'Simplex';
# arg algorithm = 'Random';
# arg time_limit = 10;
  arg total_runs = 1;
  }

  def input_params
  {
  param N[] = [10000];
  param alpha = 1;
  param beta = 1;

  decl in static double A[N][N+20] = 0;
  decl in static double B[N][N+20] = 0;
  decl in static double x[N] = 0;
  decl in static double u1[N] = 0;
  decl in static double u2[N] = 0;
  decl in static double v1[N] = 0;
  decl in static double v2[N] = 0;
  decl out static double w[N] = 0;
  decl in static double y[N] = 0;
  decl in static double z[N] = 0;
  }
) @*/

/**-- (Generated by ancc)
Input problem sizes:
  N = 10000 
  alpha = 1 
  beta = 1 
Best performance parameters:
  IVEC1 = True 
  IVEC2 = True 
  PAR = False 
  PERM1 = [1, 0] 
  PERM2 = ['i', 'j'] 
  SCREP1 = True 
  SCREP2 = True 
  T1_1 = 1 
  T1_2 = 1 
  T2_1 = 1 
  T2_2 = 1 
  U1 = 1 
  U2 = 4 
  U3 = 4 
  U4 = 1 
  U5 = 4 
--**/

   

register int i,j,k,t;  
register int it,jt,kt,tt;  
register int c1t, c2t, c3t, c4t, c5t, c6t, c7t, c8t, c9t, c10t, c11t, c12t;  
register int newlb_c1, newlb_c2, newlb_c3, newlb_c4, newlb_c5, newlb_c6,  
  newlb_c7, newlb_c8, newlb_c9, newlb_c10, newlb_c11, newlb_c12;  
register int newub_c1, newub_c2, newub_c3, newub_c4, newub_c5, newub_c6,  
  newub_c7, newub_c8, newub_c9, newub_c10, newub_c11, newub_c12;  

/*@ begin PolySyn(  
  parallel = PAR; 
  tiles = [T1_1,T1_2,T2_1,T2_2];
  permut = PERM1;
  unroll_factors = [U1,U2]; 
  scalar_replace = SCREP1; 
  vectorize = IVEC1; 
 
  profiling_code = 'gemver_profiling.c'; 
  compile_cmd = 'gcc'; 
  compile_opts = '-lm'; 
  ) @*/ 



#define ceild(n,d)  ceil(((double)(n))/((double)(d)))
#define floord(n,d) floor(((double)(n))/((double)(d)))
#define max(x,y)    ((x) > (y)? (x) : (y))
#define min(x,y)    ((x) < (y)? (x) : (y))

		
	int c1, c2, c3, c4, c5, c6, c7;

/* Generated from PLuTo-produced CLooG file by CLooG v0.14.1 64 bits in 0.01s. */
/*@ begin Loop(
transform Composite(
permut = [['c2', 'c1']],
  regtile = (['c1', 'c2'],[1, 4]),
  scalarreplace = (True, 'double'),
  vector = (True, ['ivdep','vector always']))

for (c1=0;c1<=N-1;c1++) {
  for (c2=0;c2<=N-1;c2++) {
    B[c2][c1]=u2[c2]*v2[c1]+u1[c2]*v1[c1]+A[c2][c1] ;
    x[c1]=beta*B[c2][c1]*y[c2]+x[c1] ;
  }
}

) @*/{
  for (c2t=0; c2t<=N-4; c2t=c2t+4) {
    register int cbv_1;
    cbv_1=N-1;
#pragma ivdep
#pragma vector always
    for (c1=0; c1<=cbv_1; c1++ ) {
      double scv_1, scv_2, scv_3, scv_4, scv_5, scv_6, scv_7;
      scv_1=B[c2t][c1];
      scv_2=B[(c2t+3)][c1];
      scv_3=v2[c1];
      scv_4=B[(c2t+1)][c1];
      scv_5=B[(c2t+2)][c1];
      scv_6=x[c1];
      scv_7=v1[c1];
      scv_1=u2[c2t]*scv_3+u1[c2t]*scv_7+A[c2t][c1];
      scv_4=u2[(c2t+1)]*scv_3+u1[(c2t+1)]*scv_7+A[(c2t+1)][c1];
      scv_5=u2[(c2t+2)]*scv_3+u1[(c2t+2)]*scv_7+A[(c2t+2)][c1];
      scv_2=u2[(c2t+3)]*scv_3+u1[(c2t+3)]*scv_7+A[(c2t+3)][c1];
      scv_6=beta*scv_1*y[c2t]+scv_6;
      scv_6=beta*scv_4*y[(c2t+1)]+scv_6;
      scv_6=beta*scv_5*y[(c2t+2)]+scv_6;
      scv_6=beta*scv_2*y[(c2t+3)]+scv_6;
      B[c2t][c1]=scv_1;
      B[(c2t+3)][c1]=scv_2;
      B[(c2t+1)][c1]=scv_4;
      B[(c2t+2)][c1]=scv_5;
      x[c1]=scv_6;
    }
  }
  for (c2=c2t; c2<=N-1; c2=c2+1) {
    register int cbv_2;
    cbv_2=N-1;
#pragma ivdep
#pragma vector always
    for (c1=0; c1<=cbv_2; c1++ ) {
      double scv_8, scv_9;
      scv_8=x[c1];
      scv_9=B[c2][c1];
      scv_9=u2[c2]*v2[c1]+u1[c2]*v1[c1]+A[c2][c1];
      scv_8=beta*scv_9*y[c2]+scv_8;
      x[c1]=scv_8;
      B[c2][c1]=scv_9;
    }
  }
}
/*@ end @*/

/* End of CLooG code */

/*@ end @*/

/*@ begin Loop(
 transform Composite(
  regtile = [['i'],[U5]],
  scalarreplace = (SCREP2, 'double', 'itv_'),
  vector = (IVEC2, ['ivdep','vector always'])
 )
for (i=0; i<=N-1; i++)
  x[i] = x[i] + z[i];
 ) @*/{
  register int cbv_1;
  cbv_1=N-4;
#pragma ivdep
#pragma vector always
  for (it=0; it<=cbv_1; it=it+4) {
    double itv_1, itv_2, itv_3, itv_4;
    itv_1=x[(it+3)];
    itv_2=x[(it+1)];
    itv_3=x[(it+2)];
    itv_4=x[it];
    itv_4=itv_4+z[it];
    itv_2=itv_2+z[(it+1)];
    itv_3=itv_3+z[(it+2)];
    itv_1=itv_1+z[(it+3)];
    x[(it+3)]=itv_1;
    x[(it+1)]=itv_2;
    x[(it+2)]=itv_3;
    x[it]=itv_4;
  }
  register int cbv_2;
  cbv_2=N-1;
#pragma ivdep
#pragma vector always
  for (i=it; i<=cbv_2; i=i+1) {
    double itv_5;
    itv_5=x[i];
    itv_5=itv_5+z[i];
    x[i]=itv_5;
  }
}
/*@ end @*/

/*@ begin Loop(
 transform Composite(
  permut = [PERM2],
  regtile = [['i','j'],[U3,U4]],
  scalarreplace = (SCREP2, 'double', 'itv_'),
  vector = (IVEC2, ['ivdep','vector always'])
 )
 for (i=0; i<=N-1; i++)
   for (j=0; j<=N-1; j++)
     w[i] = w[i] + alpha* B[i][j]*x[j];
 ) @*/{
  for (it=0; it<=N-4; it=it+4) {
    double itv_2, itv_3, itv_4, itv_5;
    itv_2=w[(it+2)];
    itv_3=w[(it+1)];
    itv_4=w[(it+3)];
    itv_5=w[it];
    register int cbv_1;
    cbv_1=N-1;
#pragma ivdep
#pragma vector always
    for (j=0; j<=cbv_1; j++ ) {
      double itv_1;
      itv_1=x[j];
      itv_5=itv_5+alpha*B[it][j]*itv_1;
      itv_3=itv_3+alpha*B[(it+1)][j]*itv_1;
      itv_2=itv_2+alpha*B[(it+2)][j]*itv_1;
      itv_4=itv_4+alpha*B[(it+3)][j]*itv_1;
    }
    w[(it+2)]=itv_2;
    w[(it+1)]=itv_3;
    w[(it+3)]=itv_4;
    w[it]=itv_5;
  }
  for (i=it; i<=N-1; i=i+1) {
    double itv_6;
    itv_6=w[i];
    register int cbv_2;
    cbv_2=N-1;
#pragma ivdep
#pragma vector always
    for (j=0; j<=cbv_2; j++ ) {
      itv_6=itv_6+alpha*B[i][j]*x[j];
    }
    w[i]=itv_6;
  }
}
/*@ end @*/



/*@ end @*/


    annot_t_end = rtclock();
    annot_t_total += annot_t_end - annot_t_start;
  }
  
  annot_t_total = annot_t_total / REPS;
  printf("%f\n", annot_t_total);

  return ((int) w[0]); 

}
                                    
